Exceptional population genomic homogeneity in the black brittle star Ophiocomina nigra (Ophiuroidea, Echinodermata) along the Atlantic-Mediterranean coast

The Atlantic-Mediterranean marine transition is characterised by strong oceanographic barriers and steep environmental gradients that generally result in connectivity breaks between populations from both basins and may lead to local adaptation. Here, we performed a population genomic study of the black brittle star, Ophiocomina nigra, covering most of its distribution range along the Atlantic-Mediterranean region. Interestingly, O. nigra is extremely variable in its coloration, with individuals ranging from black to yellow-orange, and different colour morphs inhabiting different depths and habitats. In this work, we used a fragment of the mitochondrial COI gene and 2,374 genome-wide ddRADseq-derived SNPs to explore: (a) whether the different colour morphs of O. nigra represent different evolutionary units; (b) the disruptive effects of major oceanographic fronts on its population structure; and (c) genomic signals of local adaptation to divergent environments. Our results revealed exceptional population homogeneity, barely affected by oceanographic fronts, with no signals of local adaptation nor genetic differentiation between colour morphs. This remarkable panmixia likely results from a long pelagic larval duration, a large effective population size and recent demographic expansions. Our study unveils an extraordinary phenotypic plasticity in O. nigra, opening further research questions on the ecological and molecular mechanisms underpinning coloration in Ophiuroidea.

The Atlantic-Mediterranean biogeographic marine region encompasses a wide range of subtropical, temperate and even subarctic marine areas that also includes several environmental transition zones and oceanographic fronts [1][2][3][4] . One of these relevant oceanographic transitions is the English Channel, which brings the southern Atlantic water mass into contact with the colder and shallower North Sea 5,6 . In the southern European area, the Strait of Gibraltar and the Almeria-Oran front define two oceanographic breaks, due to the strong and well-defined marine circulation system across this area, that separate the Atlantic and the Mediterranean basins 6,7 . These two barriers delimit the Alboran Sea, a mixing-water zone between the warmer and saltier Mediterranean Sea and the lower-salinity and colder Atlantic Ocean 8 . Depending on the species, these marine barriers have different effects on the connectivity and gene flow patterns between Atlantic and Mediterranean populations 6,7,[9][10][11][12][13][14] , thus providing different disruptive effects of these two oceanographic barriers 3,9,10,12,13,[15][16][17][18] . Genetic diversity patterns highlighted that large-and short-scale oceanographic circulation interplays with each species life history traits and its evolutionary history to determine its divergence level 3,6,19,20 . Therefore, www.nature.com/scientificreports/ From these 192 individuals, we sub-sampled 144 organisms that were used for Single Nucleotide Polymorphisms (SNP) isolation using double digest restriction-site associated DNA sequencing (ddRADseq). After filtering the Illumina reads, we obtained a total of ca. 272 M reads, and retained 109 individuals with an average of ca. 1.9 M reads per individual. The final dataset included 2,374 neutral SNPs and 109 individuals, containing 30.46% of missing data and an average depth of 9.6 reads per SNP per sample.
We detected a total of 124 outlier SNPs using Arlequin, as potential candidates under positive selection. However, no SNP was identified as under selection using Bayescan. Therefore, since we did not find common SNPs between the two methods, all SNPs were considered as neutral for subsequent population genomic analyses.   Table 2. Population genomic statistics for the SNP dataset are shown in Table 3. Allele number and allele effective number were similar across the sampling stations, ranging from 1.640 in Kristineberg (KRI) to 1.956 in Tarifa (TAR) for the number of alleles, and from 1.337 in Kristineberg (KRI) to 1.376 in Tarifa (TAR) for the effective number of alleles. Genetic diversity, measured as observed heterozygosity (Ho) ranged from 0.203 in Es Caials (CAI) to 0.233 in Cabo Negro (CNE), with a value of 0.222 for the entire dataset, and expected heterozygosity (He), ranged from 0.236 in Es Caials (CAI) to 0.253 in Cabo Negro (CNE), with a value of 0.246 for the whole dataset. Observed heterozygosity was lower than the expected value in all sites and for the whole dataset. Inbreeding coefficients (F IS ) presented low values in all sampling sites, ranging from 0.079 in Cabo Negro (CNE) to 0.138 in Es Caials (CAI), with a total value of 0.099. All populations were deviated from the Hardy-Weinberg equilibrium as demonstrated by their significant values.
Population connectivity and structure. COI results. The COI haplotype network is shown in Fig. 2 (haplotype frequencies per population available in Mendeley Data at https:// doi. org/ 10. 17632/ 5m6v5 sm5f6.1). It reveals deep diverse bush-like genealogies, with twenty shared haplotypes and a large number of singletons. The shared haplotypes suggest population homogeneity and high population connectivity, with haplotypes shared among individuals from the three different regions studied (Fig. 2). Interestingly, no structure was either detected between colour morphs, pointing to phenotypic plasticity as the mechanism explaining the environmental segregation of colour morphs in O. nigra. In agreement with these results, the pairwise ϕ ST table (Supplementary Material SM1) showed low and non-significant values for all pairwise comparisons, ranging from 0.001 to 0.032, revealing low differentiation between sampling stations across our study area. Analyses of molecular variance (AMOVAs) showed that there was no significant genetic differentiation by region or colour morph in the COI dataset (Supplementary material SM2), in agreement with the previously exposed results showing the high genetic homogeneity of O. nigra among sampling stations and regions.
ddRADseq-derived SNP results. Population structure results for the ddRADseq-derived SNP dataset also showed high levels of homogeneity across regions. Pairwise F ST table (Supplementary Material SM1) showed low to moderate values of differentiation, ranging from 0 to 0.168. Only two comparisons were significantly different from zero (p-value after correction < 0.05: Cabo Negro (CNE) -Tarifa (TAR) and Cabo Negro (CNE) -La Herradura (HER), despite the short geographic distance among these three locations). AMOVA results for  56 , showed that one cluster was again the optimal number to describe our data (Supplementary Material SM6). However, a subtle substructure emerged when sampling site information was used as a prior for the DAPC (Fig. 3B), with some differentiation between the three regions. The complete Table 3. Population genetic diversity indices for the ddRADseq-derived SNP dataset including the number of individuals retained after filtering (Nf), number of alleles (An), effective number of alleles (A_Eff n), observed heterozygosity (Ho), expected heterozygosity (He), and inbreeding index (F IS ) with the Hardy-Weinberg Equilibrium p-value (**p-value < 0.01). We also included the three neutrality tests performed: Tajima's D, Fu's Fs and R2. Data are shown by locality, colour morph and region. The Mismatch distribution indicated that the COI data did not fit an unimodal expansion model, as shown by the significant values of the r statistic (r = 0.0027, p = 0.001), but the bimodal distribution on the observed pairwise differences may correspond to two consecutive expansion events (Supplementary Material SM7). The coalescent-based Bayesian Skyline Plot (BSP) for the COI dataset supported the existence of two population growth events, the most ancient one occurring between 50,000-60,000 years ago, and the most recent starting at 10,000 years ago until the present days (Fig. 4A).
Neutrality tests for the ddRAD data are shown by locality, colour morph and region (Table 3)  www.nature.com/scientificreports/ Stairway Plot 2 57 results using the folded site frequency spectrum derived from the SNP dataset revealed the particular timing when the most recent population expansion occurred (Fig. 4B). This recent expansion started at around 50,000 years ago, in agreement with the BSP results, and presented a very steep population growth between 20,000 and 15,000 years ago.

Discussion
The study here presented is the first exploring the population genomic structure, phylogeography, and demographic history patterns of the brittle star Ophiocomina nigra. To date, genome-wide scans for population genomics of ophiuroids only focused on three Antarctic species 43,44 . Although O. nigra has been deemed as uncommon in the Mediterranean Sea 37,52 , we report here populations in our Mediterranean sampling sites. We suggest that this is not due to a recent expansion of this species into the Mediterranean Sea in the last 80 years, as we did not find strong genetic signals of founder effects. We consider that the Mediterranean populations were overlooked in the past.
Our population genomic results, based on 2,374 ddRADseq-derived SNPs and a fragment of the mitochondrial COI marker, revealed population homogeneity across the whole study area, from the cold waters of the Northeastern Atlantic Ocean to the Mediterranean Sea. Our results demonstrated that the two different colour morphs present in O. nigra correspond to a single evolutionary unit, as no genetic divergence was observed between them. These results contrast with the high number of cryptic speciation events detected in other brittle stars in the same geographical area, such as Ophiothrix spp 18,31 , Ophioderma longicauda and Amphipholis squamata species complexes 30,34 . Hence, our results point to a large phenotypic plasticity as the driver of the morphological differences found between O. nigra colour morphs. As previously suggested, colour variation in O. nigra may depend on the specific environment where an organism lives 52 . In general, individuals inhabiting shallower depths tend to exhibit darker colours, whereas in deeper habitats, the individuals have lighter tones 52 . These differences in pigmentation have been associated with the higher melanin production that individuals from shallower habitats present as a form of sun protection 52,58 . Another factor that has been suggested to affect body coloration in O. nigra is the substrate characteristics 52 . However, we cannot test this hypothesis because no details on the substrate type were collected during our sampling.
Despite this morphological variation within the species, and the large geographical range explored, we found general population homogeneity. We detected high levels of genetic diversity for both the mitochondrial COI fragment and the nuclear SNPs, but no population divergence (see Figs. 1, 2 and 3). The major oceanographic fronts included in our study area, the English Channel, the Gibraltar Strait, and the Almería-Oran front had a minor disruptive effect on the population structure of O. nigra, a result that contrasts with those observed in most coastal benthic species (e.g. 3,7,18 ), including other brittle stars 18 . Such a pattern of panmixia, or low genetic divergence, between the Mediterranean and Atlantic basins, or within the Mediterranean sub-basins is relatively common in echinoderms with large dispersal potential 15,21,59 and in those with a predominance of asexual reproduction 9,17 . Nevertheless, the absence of divergence found in O. nigra between the Northeast Atlantic and the Mediterranean is an exceptional pattern even for species with large dispersal potential 12,14,15,21,59 . Past oceanographic processes across the Atlantic-Mediterranean transition 6 , and the current disruptive effect of major oceanographic circulation across the Gibraltar Strait and the Almería-Oran front, have generated a large variety of genomic patterns in marine invertebrates, but most of them are characterised by a major divergence between these two main basins 12,14,16 , although exceptions also exist 31 . Interestingly, within the Ophiothrix species complex, Ophiothrix sp. II displayed a lack of genetic divergence similar to that obtained in O. nigra, with genetic www.nature.com/scientificreports/ homogeneity along the Atlantic-Mediterranean transition 18,31 . However, for Ophiothrix sp. II the sampling area was restricted to the Iberian Peninsula and only two mitochondrial markers were studied 18 . The lack of significant differentiation between distant geographical areas in O. nigra may result from the combination of different biological features, demographic and historical events, including a large dispersal potential of the larva, large effective population size, and past demographic expansions 14,18,31 . The planktotrophic larva of O. nigra remains in the water column for several months avoiding settlement near adults 51 , and hence can be dispersed hundreds of kilometres by marine currents 47 . A long dispersal potential during early life stages, which interacts with oceanographic circulation, significantly boosts population connectivity 60 . Additionally, the low-density and regularly dispersed patches formed by O. nigra across large spatial areas 47 may favour connectivity along distant sites.
The population homogeneity found in O. nigra could be also the result of the recent demographic expansions that we identified from both COI and SNPs datasets. Our demographic history analyses identified a recent effective population expansion consisting in two growth periods, likely related to expansions during Pleistocene interglacial periods (e.g. 12,18 ). The ancient growth period was recovered from both the COI BSP and the ddRADseq Stairway Plot 2 results, and was dated between 50,000 and 60,000 years ago, similar to that dated for Ophiothrix spp 18 and some other Mediterranean species 61 . The most recent one was also identified by both datasets, and pointed to a tenfold increase of O. nigra populations between 10,000 and 20,000 years ago (see Fig. 4). This coincides with the deglaciation period that followed the Last Glacial Maximum (LGM) 62 , and is the most common post-LGM population expansion pattern found in other marine invertebrates (see an example in 61 ). This recent demographic expansion could have partially erased previous population differentiation in O. nigra, favouring the current genetic homogeneity found in both datasets. Additionally, a large effective population size reduces genetic drift effects 63 , reducing genetic differentiation between populations and increasing genetic homogeneity.
Despite the general pattern of population homogeneity found in O. nigra, the presence of sub-regional private haplotypes in the COI marker suggested a certain degree of genetic substructuring in major geographical areas. In agreement with these results, migration analyses based on ddRADseq-derived SNPs showed that gene flow is not symmetric and homogeneous over our study area. The migration network (Fig. 3C) showed that Kristineberg (KRI) was isolated from all other sites, whereas Roscoff (ROS) and the NW Mediterranean sites were all connected to two sites from the Alboran Sea, Tarifa (TAR) and La Herradura (HER). This asymmetric migration pattern reveals a disruptive effect of the Gibraltar Strait or the Almería-Oran front on gene flow. Therefore, the role of the Alboran Sea as a transitional area connecting these two major basins, the Mediterranean Sea and the Atlantic Ocean 3,9,15,22 , is also found in O. nigra.
Finally, against our initial expectations, we did not find genomic signals of selection and/or local adaptation in O. nigra across our studied area, which covers important temperature and salinity clines. This result in O. nigra contrasts with those from other echinoderms that displayed signals of selection and local adaptation to temperature and salinity 15,21 , even under high levels of gene flow. However, due to the small proportion of the genome that is covered by RADseq approaches, strong selection associated with environmental variables might be affecting other genomic regions that remained uncovered in this study 64 .
Overall, our study shows an exceptionally high genetic homogeneity in O. nigra, barely disturbed by the disruptive effects of the English Channel, the Gibraltar Strait and the Almería-Oran front. This genetic homogeneity likely results from a combination of a long pelagic larval duration, a large effective population size, and recent demographic expansions. Our study also demonstrates that the different colour morphs present in O. nigra are the result of high phenotypic plasticity of the species, opening new questions for further research investigating the ecological and molecular mechanisms behind this intriguing plasticity.

Methods
Sampling, DNA extraction and sequencing. A total of 192 specimens of Ophiocomina nigra were sampled across the natural distribution range of the species, from 10 different localities throughout the Atlantic-Mediterranean region between 2014 and 2017 (see details in Fig. 1 and Table 1). Specimens were collected between 5-20 m depth by SCUBA diving for most localities, except for those from Blanes (BLA, NW Mediterranean), which were collected by local fishermen trawling between 100-150 m depth. All samples were preserved in absolute ethanol at room temperature until further processing.
Among all samples, we identified the two colour morphs previously described for this species and likely associated with depth (see Table 3 in 19 ). See collection details in Table 1. www.nature.com/scientificreports/ ddRAD library preparation, sequencing and filtering. A subset of 144 specimens from eight representative O. nigra populations, was used to generate double digest restriction site associated DNA (ddRAD) libraries following 67 . High molecular weight DNA was re-extracted from preserved samples following a modified cetyltrimethylammoniumbromide (CTAB) protocol 68 . Briefly, tube feet tissue samples were digested overnight at 57 °C in CTAB buffer (2% w/v CTAB, 1.4 M NaCl, 20 mM EDTA, 100 mM Tris-HCl, pH 8), 2-mercaptoethanol and Proteinase K. DNA was subsequently isolated with 24:1 chloroform:isoamyl alcohol and precipitated with ice-cold ethanol. DNA integrity, purity, and quantity were checked in agarose electrophoresis gels, Nanodrop (Thermo Fisher), and using a Qubit DNA HS assay (Life Technologies), respectively. For ddRAD library preparation, we double digested 510 ng of DNA per sample using high-fidelity restriction enzymes EcoRI and MseI (New England Biolabs). Fragments per sample were purified with Agencourt AMPure beads (Beckham Coulter) before ligation with barcoded Illumina adapters. Samples with unique adapters were pooled, and each pool was size-selected (200 to 400 bp) using a Blue Pippin (Sage Science). Illumina multiplexing indices were added to each library during PCR amplification with a Phusion high-fidelity DNA polymerase kit (New England Biolabs), using a common forward primer (PCR1) and a reverse library-specific primer (indexed PCR2). Six libraries, with 24 individuals per library, were paired-end sequenced (150 bp) on an Illumina NovaSeq 6000 at the Novogene Core Facility, Cambridgeshire, UK. Quality filtering, locus assembly and SNP calling were conducted using the Stacks pipeline v. 2.59 69 . RADtags (DNA fragments with the two appropriate restriction enzyme cut sites that were selected, amplified, and sequenced) were processed using process_radtags, where raw reads were quality-trimmed to remove low quality reads, reads with uncalled bases, and reads without a complete barcode or restriction cut site. After these filtering steps, we retained a total of 272,024,358 reads from the initial 307,074,382 raw reads, with an average of 1,889,058 reads per sample. Then, the core Stacks pipeline was run with optimised parameters following 70 and 71 . Optimal parameter values were: m = 3, M = 4, and n = 4. The Stacks population module was used to conduct a first filtering of the data, retaining SNPs present in at least 50% of the individuals (r = 0.5) and with a minimum allele frequency higher than 0.05 (-min-maf 0.05). We also used the "-write-single-snp" option in order to retain only the first SNP from each RAD-tag to prevent the inclusion of physically linked loci in the dataset. We then filtered out individuals with less than 20% of the loci using the adegenet R package 72 . Then, we used the perl script "filter_ hwe_by_pop.pl" from the dDocent pipeline 73 with default parameters (-h 0.001, -c 0.25) to filter SNPs deviating from Hardy-Weinberg equilibrium within sampling stations.
The identification of putative SNPs under selection was performed using Arlequin v. 3.5 74 and Bayescan v. 2.1 75 . Arlequin was run with default parameters using 100,000 coalescent simulations and 100 demes. Benjamini-Yekutieli FDR corrections were applied to 0.01 p-value. Bayescan was run with default parameters, using the -snp flag and a false discovery rate of 0.05 (FDR = 0.05). Only SNPs identified by both methods were considered as being under selection and removed from the neutral SNP dataset used for subsequent population genomic analyses.
Population structure and genetic diversity analyses. Mitochondrial COI gene fragment. We calculated diversity indices for the COI data, including haplotype diversity (Hd), nucleotide diversity (π), the number of haplotypes for the different populations (H), and number of private haplotypes per sampling site, colour morph and region using DnaSP v. 5.10 76 . The software PopART 77 was used to construct an un-rooted haplotype network using the TCS statistical parsimony network approach 78 .
Pairwise F ST indices were estimated between sampling stations using Arlequin with 90,600 permutations to determine their associated p-values. We used a Benjamini-Yekutieli false discovery rate correction for multiple comparisons 79 and an overall corrected 0.05 α-level. In order to test whether the groupings of samples by colour morph or by region explained a significant part of the total genetic variation we performed two analyses of molecular variance (AMOVA) in Arlequin.
Neutral ddRAD data. Indices of genetic diversity including mean number of alleles, expected heterozygosity (He), observed heterozygosity (Ho), inbreeding index (F IS ), and the Hardy-Weinberg Equilibrium were estimated by sampling site, colour morph and region using GenoDive v. 3.05 80 .
Pairwise F ST indices and AMOVAs were also estimated for the ddRADseq-derived SNP dataset using Arlequin as detailed above. Population structure was assessed using STRU CTU RE v. 2.3.4 55 and the discriminant analysis of principal components (DAPC) 56 as implemented in the adegenet R package 72 . STRU CTU RE was run using the admixture model for 500,000 MCMC iterations, with a burn-in of 50,000 iterations, setting the number of clusters (K) from 1 to 9, with 10 replicates per K. We used STRU CTU RE HARVESTER 81 to inspect the most likely number of K, and CLUMPAK 82 to average individual's membership coefficients across replicates and graphically represent STRU CTU RE results. The DAPC was performed estimating the number of clusters with the find.clusters function in adegenet, using the optim.a.score function to optimise the number of retained PCs. Additionally, in order to identify gene flow patterns, a relative migration network analysis was conducted with the divMigrate function in the diveRsity R package 83,84 using Nei's G ST method.
Demographic history analyses. In order to unveil the evolutionary history of Ophiocomina nigra and infer how it might have affected the current genetic diversity distribution and population structure, we ran demographic history analyses for both COI and ddRAD data.
For both datasets, Tajima's D and Ramos-Onsins and Rozas' R2 neutrality tests were calculated for each sampling station, colour morph and region, and for the whole datasets using DnaSP v. 6.12 85  www.nature.com/scientificreports/ sampling site, colour morph and region, using DnaSP. We used the fasta file generated by the Stacks population module with the -fasta-samples option as the ddRAD data input file for DnaSP. For the COI data, we calculated the pairwise mismatch distribution using DnaSP and performed a Bayesian Skyline Plot (BSP) using BEAUTi and BEAST v. 1.10.4 86 . For the pairwise mismatch distribution, the effective population size parameter (θ) was set to follow the Rogers & Harpending model 87 . In order to determine if the empirical data fitted an expansion model we calculated the raggedness (r) index 88 using 10,000 permutations. Our priors for the BSP model included the implementation of the substitution model previously defined by jModelTest v. 2.1.10 89 (HKY + I + G), a strict clock model, and the constant skyline model. As no molecular clock has been calibrated for ophiuroids, a mutation rate of 1.01% per nucleotide per million years was used for the COI, following other studies in echinoderms and the most recent estimations for this group 90 . Analyses were run for 500 million generations, with a sampling frequency of 50,000 generations, and a burning of 25% of the MCMC. The software Tracer v.1.7.2 91 was used for assessing stationarity of the MCMC, effective sample sizes (ESSs > 200) and generating the evolution of effective population size under the skyline plot model, expressed as NeT (T = generation time) over time.